C
C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C LINK-MT3DMS (LMT) PACKAGE V7 FOR MODFLOW-2005
C Modified from LMT V6 for MODFLOW-2000 as documented in:
C     Zheng, C., M.C. Hill, and P.A. Hsieh, 2001,
C         MODFLOW-2000, the U.S. Geological Survey modular ground-water
C         model--User guide to the LMT6 Package, the linkage with
C         MT3DMS for multispecies mass transport modeling:
C         U.S. Geological Survey Open-File Report 01-82
C
C Revision History: 
C     Version 7.0: 08-08-2008 cz
C     Version 7.0: 08-15-2009 swm: added LMTMODULE to support LGR
C     Version 7.0: 02-12-2010 swm: rolled in include file
C ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C
C
      MODULE LMTMODULE
         INTEGER, SAVE, POINTER ::ISSMT3D,IUMT3D,ILMTFMT
        TYPE LMTTYPE
         INTEGER, POINTER ::ISSMT3D,IUMT3D,ILMTFMT
        END TYPE
        TYPE(LMTTYPE), SAVE  ::LMTDAT(10)
      END MODULE LMTMODULE
C
      SUBROUTINE LMT7BAS7AR(INUNIT,CUNIT,IGRID)
C **********************************************************************
C OPEN AND READ THE INPUT FILE FOR THE LINK-MT3DMS PACKAGE VERSION 7.
C CHECK KEY FLOW MODEL INFORMATION AND SAVE IT IN THE HEADER OF
C THE MODFLOW-MT3DMS LINK FILE FOR USE IN MT3DMS TRANSPORT SIMULATION.
C NOTE THE 'STANDARD' HEADER OPTION IS NO LONGER SUPPORTED. INSTEAD,
C THE 'EXTENDED' HEADER OPTION IS THE DEFAULT. THE RESULTING LINK FILE 
C IS ONLY COMPATIBLE WITH MT3DMS VERSION [4.00] OR LATER.
C **********************************************************************
C last modified: 08-08-2008
C last modified: 10-21-2010 swm: added MTMNW1 & MTMNW2
C      
      USE GLOBAL,   ONLY:NCOL,NROW,NLAY,NPER,NODES,NIUNIT,IUNIT,
     &                   ISSFLG,IBOUND,IOUT
      USE LMTMODULE,ONLY:ISSMT3D,IUMT3D,ILMTFMT 

C--USE FILE SPECIFICATION of MODFLOW-2005
      INCLUDE 'openspec.inc'
      LOGICAL       LOP
      CHARACTER*4   CUNIT(NIUNIT)
      CHARACTER*200 LINE,FNAME,NME
      CHARACTER*8   OUTPUT_FILE_HEADER
      CHARACTER*11  OUTPUT_FILE_FORMAT      
      DATA          INLMT,MTBCF,MTLPF,MTHUF,MTWEL,MTDRN,MTRCH,MTEVT,
     &              MTRIV,MTSTR,MTGHB,MTRES,MTFHB,MTDRT,MTETS,MTSUB,
     &              MTIBS,MTLAK,MTMNW,MTSWT,MTSFR,MTUZF
     &             /22*0/
C     -----------------------------------------------------------------    
      ALLOCATE(ISSMT3D,IUMT3D,ILMTFMT)
C
C--SET POINTERS FOR THE CURRENT GRID 
cswm: already set in main (GWF2BAS7OC)      CALL SGWF2BAS7PNT(IGRID)     
C
C--CHECK for OPTIONS/PACKAGES USED IN CURRENT SIMULATION
      IUMT3D=0
      DO IU=1,NIUNIT
        IF(CUNIT(IU).EQ.'LMT6') THEN
          INLMT=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'BCF6') THEN
          MTBCF=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'LPF ') THEN
          MTLPF=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'HUF2') THEN
          MTHUF=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'WEL ') THEN
          MTWEL=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'DRN ') THEN
          MTDRN=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'RCH ') THEN
          MTRCH=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'EVT ') THEN
          MTEVT=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'RIV ') THEN
          MTRIV=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'STR ') THEN
          MTSTR=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'GHB ') THEN
          MTGHB=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'RES ') THEN
          MTRES=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'FHB ') THEN
          MTFHB=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'DRT ') THEN
          MTDRT=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'ETS ') THEN
          MTETS=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'SUB ') THEN
          MTSUB=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'IBS ') THEN
          MTIBS=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'LAK ') THEN
          MTLAK=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'MNW1') THEN
          MTMNW1=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'MNW2') THEN
!swm: store seperate to not get clobbered by MNW1
          MTMNW2=IUNIT(IU)   
        ELSEIF(CUNIT(IU).EQ.'SWT ') THEN
          MTSWT=IUNIT(IU)        
        ELSEIF(CUNIT(IU).EQ.'SFR ') THEN
          MTSFR=IUNIT(IU)
        ELSEIF(CUNIT(IU).EQ.'UZF ') THEN
          MTUZF=IUNIT(IU)
        ENDIF
      ENDDO            
!swm: SET MTMNW IF EITHER MNW1 OR MNW2 IS ACTIVE
      IF(MTMNW1.NE.0) MTMNW=MTMNW1
      IF(MTMNW2.NE.0) MTMNW=MTMNW2
C
C--IF LMT7 PACKAGE IS NOT ACTIVATED, SKIP TO END AND RETURN
      IF(INLMT.EQ.0) GOTO 9999
C
C--ASSIGN DEFAULTS TO LMT INPUT VARIABLES AND OUTPUT FILE NAME
      IUMT3D=333
      OUTPUT_FILE_HEADER='EXTENDED'
      ILMTHEAD=1
      OUTPUT_FILE_FORMAT='UNFORMATTED'
      ILMTFMT=0     
      INQUIRE(UNIT=INLMT,NAME=NME,OPENED=LOP)
      IFLEN=INDEX(NME,' ')-1
      DO NC=IFLEN,2,-1
        IF(NME(NC:NC).EQ.'.') THEN      
          FNAME=NME(1:NC-1)//'.FTL'
          GO TO 10
        ENDIF
      ENDDO    
      FNAME=NME(1:IFLEN)//'.FTL'     
C
C--READ ONE LINE OF LMT PACKAGE INPUT FILE
   10 READ(INLMT,'(A)',END=1000) LINE
      IF(LINE.EQ.' ') GOTO 10
      IF(LINE(1:1).EQ.'#') GOTO 10
C
C--DECODE THE INPUT RECORD
      LLOC=1
      CALL URWORD(LINE,LLOC,ITYP1,ITYP2,1,N,R,IOUT,INLMT)
C
C--CHECK FOR "OUTPUT_FILE_NAME" KEYWORD AND GET FILE NAME
      IF(LINE(ITYP1:ITYP2).EQ.'OUTPUT_FILE_NAME') THEN
        CALL URWORD(LINE,LLOC,INAM1,INAM2,0,N,R,IOUT,INLMT)
        IFLEN=INAM2-INAM1+1
        IF(LINE(INAM1:INAM2).EQ.' ') THEN
        ELSE
          FNAME=LINE(INAM1:INAM2)
        ENDIF
C
C--CHECK FOR "OUTPUT_FILE_UNIT" KEYWORD AND GET UNIT NUMBER
      ELSEIF(LINE(ITYP1:ITYP2).EQ.'OUTPUT_FILE_UNIT') THEN
        CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IU,R,IOUT,INLMT)
        IF(IU.GT.0) THEN
          IUMT3D=IU
        ELSEIF(IU.LT.0) THEN
          WRITE(IOUT,11) IU
          WRITE(*,11) IU
          CALL USTOP(' ')
        ENDIF
C
C--CHECK FOR "OUTPUT_FILE_HEADER" KEYWORD AND GET INPUT VALUE
      ELSEIF(LINE(ITYP1:ITYP2).EQ.'OUTPUT_FILE_HEADER') THEN
        CALL URWORD(LINE,LLOC,ISTART,ISTOP,1,N,R,IOUT,INLMT)
        IF(LINE(ISTART:ISTOP).EQ.' '.OR.
     &     LINE(ISTART:ISTOP).EQ.'EXTENDED') THEN
          OUTPUT_FILE_HEADER='EXTENDED'
          ILMTHEAD=1
        ELSEIF(LINE(ISTART:ISTOP).EQ.'STANDARD') THEN
          WRITE(IOUT,120)
          WRITE(*,120)                   
        ELSE
          WRITE(IOUT,12) LINE(ISTART:ISTOP)
          WRITE(*,12) LINE(ISTART:ISTOP)
          CALL USTOP(' ')
        ENDIF
C
C--CHECK FOR "OUTPUT_FILE_FORMAT" KEYWORD AND GET INPUT VALUE
      ELSEIF(LINE(ITYP1:ITYP2).EQ.'OUTPUT_FILE_FORMAT') THEN
        CALL URWORD(LINE,LLOC,ISTART,ISTOP,1,N,R,IOUT,INLMT)
        IF(LINE(ISTART:ISTOP).EQ.' '.OR.
     &     LINE(ISTART:ISTOP).EQ.'UNFORMATTED') THEN
          OUTPUT_FILE_FORMAT='UNFORMATTED'
          ILMTFMT=0
        ELSEIF(LINE(ISTART:ISTOP).EQ.'FORMATTED') THEN
          OUTPUT_FILE_FORMAT='FORMATTED'
          ILMTFMT=1
        ELSE
          WRITE(IOUT,14) LINE(ISTART:ISTOP)
          WRITE(*,14) LINE(ISTART:ISTOP)
          CALL USTOP(' ')
        ENDIF
C
C--ERROR DECODING LMT INPUT KEYWORDS
      ELSE
        WRITE(IOUT,28) LINE
        WRITE(*,28) LINE
        CALL USTOP(' ')
      ENDIF
C
C--CONTINUE TO THE NEXT INPUT RECORD IN LMT FILE
      GOTO 10
C
   11 FORMAT(/1X,'ERROR READING LMT PACKAGE INPUT DATA:',
     & /1X,'INVALID OUTPUT FILE UNIT: ',I5)
   12 FORMAT(/1X,'ERROR READING LMT PACKAGE INPUT DATA:',
     & /1X,'INVALID OUTPUT_FILE_HEADER CODE: ',A)
   14 FORMAT(/1X,'ERROR READING LMT PACKAGE INPUT DATA:',
     & /1X,'INVALID OUTPUT_FILE_FORMAT SPECIFIER: ',A)
   28 FORMAT(/1X,'ERROR READING LMT PACKAGE INPUT DATA:',
     & /1X,'UNRECOGNIZED KEYWORD: ',A)
  120 FORMAT(/1X,'WARNING READING LMT PACKAGE INPUT DATA:',
     &       /1X,'[STANDARD] HEADER NO LONGER SUPPORTED; ',
     &           '[EXTENDED] HEADER USED INSTEAD.')     
C     
 1000 CONTINUE     
C
C--ENSURE A UNIQUE UNIT NUMBER FOR LINK-MT3DMS OUTPUT FILE
      IF(IUMT3D.EQ.IOUT .OR. IUMT3D.EQ.INUNIT) THEN
        WRITE(IOUT,1010) IUMT3D
        WRITE(*,1010) IUMT3D
        CALL USTOP(' ')
      ELSE
        DO IU=1,NIUNIT       
          IF(IUMT3D.EQ.IUNIT(IU)) THEN
            WRITE(IOUT,1010) IUMT3D
            WRITE(*,1010) IUMT3D
            CALL USTOP(' ')
          ENDIF
        ENDDO
      ENDIF  
 1010 FORMAT(/1X,'ERROR IN LMT PACKAGE INPT DATA:'
     &       /1X,'UNIT NUMBER GIVEN FOR FLOW-TRANSPORT LINK FILE:', 
     &        I4,' ALREADY IN USE;' 
     &       /1X,'SPECIFY A UNIQUE UNIT NUMBER.')   
C
C--OPEN THE LINK-MT3DMS OUTPUT FILE NEEDED BY MT3DMS
C--AND PRINT AN IDENTIFYING MESSAGE IN MODFLOW OUTPUT FILE  
      INQUIRE(UNIT=IUMT3D,OPENED=LOP)
      IF(LOP) THEN
        REWIND (IUMT3D)
      ELSE
        IF(ILMTFMT.EQ.0) THEN
          OPEN(IUMT3D,FILE=FNAME,FORM=FORM,ACCESS=ACCESS,
     &      ACTION=ACTION(2),STATUS='REPLACE')
        ELSEIF(ILMTFMT.EQ.1) THEN
          OPEN(IUMT3D,FILE=FNAME,FORM='FORMATTED',ACTION=ACTION(2),
     &      STATUS='REPLACE',DELIM='APOSTROPHE')
        ENDIF
      ENDIF
C
      WRITE(IOUT,30) FNAME,IUMT3D,
     &               OUTPUT_FILE_FORMAT,OUTPUT_FILE_HEADER
   30 FORMAT(//1X,'***Link-MT3DMS Package v7***',
     &        /1x,'OPENING LINK-MT3DMS OUTPUT FILE: ',A,
     &        /1X,'ON UNIT NUMBER: ',I5,
     &        /1X,'FILE TYPE: ',A,
     &        /1X,'HEADER OPTION: ',A,
     &        /1X,'***Link-MT3DMS Package v7***',/1X)
C
C--GATHER AND CHECK KEY FLOW MODEL INFORMATION
      ISSMT3D=1    !loop through all stress periods        
      DO N=1,NPER    !to check if any transient sp is used
        IF(ISSFLG(N).EQ.0) THEN
          ISSMT3D=0
          EXIT
        ENDIF
      ENDDO                  
      MTISS=ISSMT3D
      MTNPER=NPER 
C
      MTCHD=0    !loop through the entire grid to get
      DO K=1,NLAY    !total number of constant-head cells
        DO I=1,NROW
          DO J=1,NCOL
            IF(IBOUND(J,I,K).LT.0) MTCHD=MTCHD+1
          ENDDO
        ENDDO
      ENDDO
C
C--ERROR CHECKING BEFORT OUTPUT
      IF(MTEVT.GT.0.AND.MTETS.GT.0) THEN
        WRITE(IOUT,1300)
        WRITE(*,1300)
        CALL USTOP(' ')
      ENDIF    
 1300 FORMAT(/1X,'ERROR IN LMT PACKAGE INPT DATA:'
     &  /1X,'Both EVT and ETS Packages are used in flow simulation;'
     &  /1X,'Only one is allowed in the same transport simulation.')
C
C--WRITE A HEADER TO MODFLOW-MT3DMS LINK FILE
      IF(OUTPUT_FILE_HEADER.EQ.'EXTENDED') THEN        
        IF(ILMTFMT.EQ.0) THEN
          WRITE(IUMT3D) 'MT3D4.00.00',
     &     MTWEL,MTDRN,MTRCH,MTEVT,MTRIV,MTGHB,MTCHD,MTISS,MTNPER,
     &     MTSTR,MTRES,MTFHB,MTDRT,MTETS,MTSUB,MTIBS,MTLAK,MTMNW,
     &     MTSWT,MTSFR,MTUZF
        ELSEIF(ILMTFMT.EQ.1) THEN
          WRITE(IUMT3D,*) 'MT3D4.00.00',
     &     MTWEL,MTDRN,MTRCH,MTEVT,MTRIV,MTGHB,MTCHD,MTISS,MTNPER,
     &     MTSTR,MTRES,MTFHB,MTDRT,MTETS,MTSUB,MTIBS,MTLAK,MTMNW,
     &     MTSWT,MTSFR,MTUZF
        ENDIF
      ENDIF

C------SAVE POINTER DATA TO ARRARYS
      CALL SLMT7PSV(IGRID)
C
C--NORMAL RETURN
 9999 RETURN
      END
C
      SUBROUTINE LMT7BD(KKSTP,KKPER,IGRID)
C **********************************************************************
C  WRITE TERMS TO THE FLOW-TRANSPORT LINK FILE FOR USE BY MT3DMS FOR
C  TRANSPORT SIMULATIONS.  THE CODE BELOW IS COPIED FROM THE lmt7.inc 
C  FILE.  INSTEAD OF USING THE INCLUDE FILE, THE CODE IS PUT INTO THIS
C  SUBROUTINE AND CALLED FROM MAIN.
C **********************************************************************
C last modified: 02-12-2010
C     
      USE GLOBAL,ONLY:IOUT,IUNIT
      USE LMTMODULE,ONLY:ISSMT3D,IUMT3D,ILMTFMT 

C--SWM: SWAP POINTERS FOR LMT DATA TO CURRENT GRID
        CALL SLMT7PNT(IGRID)
C
C--WRITE A NOTIFICATION LINE TO MODFLOW OUTPUT FILE
        WRITE(IOUT,9876) IUMT3D,KKSTP,KKPER
 9876   FORMAT(/1X,'SAVING SATURATED THICKNESS AND FLOW TERMS ON UNIT',
     &   I5,' FOR MT3DMS',/1X,'BY THE LINK-MT3DMS PACKAGE V7',
     &   ' AT TIME STEP',I5,', STRESS PERIOD',I5/)
C
C--COLLECT AND SAVE ALL RELEVANT FLOW MODEL INFORMATION
        IF(IUNIT(1) .GT.0) 
     &   CALL LMT7BCF7(ILMTFMT,ISSMT3D,IUMT3D,KKSTP,KKPER,IGRID)
        IF(IUNIT(23).GT.0) 
     &   CALL LMT7LPF7(ILMTFMT,ISSMT3D,IUMT3D,KKSTP,KKPER,IGRID)
        IF(IUNIT(37).GT.0) 
     &   CALL LMT7HUF7(ILMTFMT,ISSMT3D,IUMT3D,
     &   KKSTP,KKPER,IUNIT(47),IGRID)
        IF(IUNIT(2) .GT.0) 
     &   CALL LMT7WEL7(ILMTFMT,IUMT3D,KKSTP,KKPER,IGRID)
        IF(IUNIT(3) .GT.0) 
     &   CALL LMT7DRN7(ILMTFMT,IUMT3D,KKSTP,KKPER,IGRID)
        IF(IUNIT(8) .GT.0) 
     &   CALL LMT7RCH7(ILMTFMT,IUMT3D,KKSTP,KKPER,IGRID)
        IF(IUNIT(5) .GT.0) 
     &   CALL LMT7EVT7(ILMTFMT,IUMT3D,KKSTP,KKPER,IGRID)
        IF(IUNIT(39).GT.0) 
     &   CALL LMT7ETS7(ILMTFMT,IUMT3D,KKSTP,KKPER,IGRID)
        IF(IUNIT(4) .GT.0) 
     &   CALL LMT7RIV7(ILMTFMT,IUMT3D,KKSTP,KKPER,IGRID)
        IF(IUNIT(7) .GT.0) 
     &   CALL LMT7GHB7(ILMTFMT,IUMT3D,KKSTP,KKPER,IGRID)
        IF(IUNIT(18).GT.0) 
     &   CALL LMT7STR7(ILMTFMT,IUMT3D,KKSTP,KKPER,IGRID)
        IF(IUNIT(17).GT.0) 
     &   CALL LMT7RES7(ILMTFMT,IUMT3D,KKSTP,KKPER,IGRID)
        IF(IUNIT(16).GT.0) 
     &   CALL LMT7FHB7(ILMTFMT,IUMT3D,KKSTP,KKPER,IGRID)
        IF(IUNIT(50).GT.0) 
     &   CALL LMT7MNW27(ILMTFMT,IUMT3D,KKSTP,KKPER,IGRID)
        IF(IUNIT(52).GT.0) 
     &   CALL LMT7MNW17(ILMTFMT,IUMT3D,KKSTP,KKPER,IGRID)
        IF(IUNIT(40).GT.0) 
     &   CALL LMT7DRT7(ILMTFMT,IUMT3D,KKSTP,KKPER,IGRID)
C
      RETURN
      END
C
      SUBROUTINE LMT7BCF7(ILMTFMT,ISSMT3D,IUMT3D,KSTP,KPER,IGRID)
C *********************************************************************
C SAVE SATURATED CELL THICKNESS; FLOW ACROSS THREE CELL INTERFACES;
C TRANSIENT FLUID-STORAGE; AND LOCATIONS AND FLOW RATES OF
C CONSTANT-HEAD CELLS FOR USE BY MT3D.  THIS SUBROUTINE IS CALLED
C ONLY IF THE 'BCF' PACKAGE IS USED IN MODFLOW.
C *********************************************************************
C Modified from Harbaugh (2005)
C last modified: 08-08-2008
C
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,ISSFLG,IBOUND,HNEW,HOLD,
     &                      BUFF,CR,CC,CV,BOTM,LBOTM
      USE GWFBASMODULE,ONLY:DELT
      USE GWFBCFMODULE,ONLY:LAYCON,SC1,SC2
      CHARACTER*16 TEXT     
      DOUBLE PRECISION HD
C
C--SET POINTERS FOR THE CURRENT GRID     
cswm: already set in GWF2BCF7BDS      CALL SGWF2BCF7PNT(IGRID)      
C      
C--GET STEADY-STATE FLAG FOR THE CURRENT STRESS PERIOD               
      ISSCURRENT=ISSFLG(KPER)
C
C--CALCULATE AND SAVE SATURATED THICKNESS
      TEXT='THKSAT'
      ZERO=0.
      ONE=1.
C
C--INITIALIZE BUFF ARRAY WITH 1.E30 FOR INACTIVE CELLS
C--OR FLAG -111 FOR ACTIVE CELLS   
      FlagInactive=1.E30
      FlagActive=-111.
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            IF(IBOUND(J,I,K).EQ.0) THEN
              BUFF(J,I,K)=FlagInactive
            ELSE
              BUFF(J,I,K)=FlagActive
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--CALCULATE SATURATED THICKNESS FOR UNCONFINED/CONVERTIBLE
C--LAYERS AND STORE IN ARRAY BUFF
      DO K=1,NLAY
        IF(LAYCON(K).EQ.0 .OR. LAYCON(K).EQ.2) CYCLE
        DO I=1,NROW
          DO J=1,NCOL
            IF(IBOUND(J,I,K).NE.0) THEN
              TMP=HNEW(J,I,K)
              BUFF(J,I,K)=TMP-BOTM(J,I,LBOTM(K))
              IF(LAYCON(K).EQ.3) THEN
                THKLAY=BOTM(J,I,LBOTM(K)-1)-BOTM(J,I,LBOTM(K))
                IF(BUFF(J,I,K).GT.THKLAY) BUFF(J,I,K)=THKLAY
              ENDIF
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--SAVE THE CONTENTS OF THE BUFFER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
        WRITE(IUMT3D) BUFF
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
        WRITE(IUMT3D,*) BUFF
      ENDIF
C
C--CALCULATE AND SAVE FLOW ACROSS RIGHT FACE
      NCM1=NCOL-1
      IF(NCM1.LT.1) GO TO 405
      TEXT='QXX'
C
C--CLEAR THE BUFFER     
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            BUFF(J,I,K)=ZERO
          ENDDO
        ENDDO
      ENDDO
C
C--FOR EACH CELL
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCM1
            IF(IBOUND(J,I,K).NE.0.AND.IBOUND(J+1,I,K).NE.0) THEN
              HDIFF=HNEW(J,I,K)-HNEW(J+1,I,K)
              BUFF(J,I,K)=HDIFF*CR(J,I,K)
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--RECORD CONTENTS OF BUFFER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
        WRITE(IUMT3D) BUFF
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
        WRITE(IUMT3D,*) BUFF
      ENDIF
C
  405 CONTINUE
C
C--CALCULATE AND SAVE FLOW ACROSS FRONT FACE
      NRM1=NROW-1
      IF(NRM1.LT.1) GO TO 505
      TEXT='QYY'
C
C--CLEAR THE BUFFER
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            BUFF(J,I,K)=ZERO
          ENDDO
        ENDDO
      ENDDO
C
C--FOR EACH CELL
      DO K=1,NLAY
        DO I=1,NRM1
          DO J=1,NCOL
            IF(IBOUND(J,I,K).NE.0.AND.IBOUND(J,I+1,K).NE.0) THEN
              HDIFF=HNEW(J,I,K)-HNEW(J,I+1,K)
              BUFF(J,I,K)=HDIFF*CC(J,I,K)
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--RECORD CONTENTS OF BUFFER.
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
        WRITE(IUMT3D) BUFF
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
        WRITE(IUMT3D,*) BUFF
      ENDIF
C
  505 CONTINUE
C
C--CALCULATE AND SAVE FLOW ACROSS FRONT FACE
      NLM1=NLAY-1
      IF(NLM1.LT.1) GO TO 700
      TEXT='QZZ'
C
C--CLEAR THE BUFFER
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            BUFF(J,I,K)=ZERO
          ENDDO
        ENDDO
      ENDDO
C
C--FOR EACH CELL CALCULATE FLOW THRU LOWER FACE & STORE IN BUFFER
      DO K=1,NLM1
        DO I=1,NROW
          DO J=1,NCOL
            IF(IBOUND(J,I,K).NE.0.AND.IBOUND(J,I,K+1).NE.0) THEN
              HD=HNEW(J,I,K+1)
              IF(LAYCON(K+1).EQ.3 .OR. LAYCON(K+1).EQ.2) THEN
                TMP=HD
                IF(TMP.LT.BOTM(J,I,LBOTM(K+1)-1))
     &           HD=BOTM(J,I,LBOTM(K+1)-1)
              ENDIF
              HDIFF=HNEW(J,I,K)-HD
              BUFF(J,I,K)=HDIFF*CV(J,I,K)
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--RECORD CONTENTS OF BUFFER.
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
        WRITE(IUMT3D) BUFF
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
        WRITE(IUMT3D,*) BUFF
      ENDIF
C
  700 CONTINUE
C
C--CALCULATE AND SAVE GROUNDWATER STORAGE IF TRANSIENT
      IF(ISSMT3D.NE.0) GO TO 705
      TEXT='STO'
C
C--INITIALIZE AND CLEAR BUFFER
      TLED=ONE/DELT
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            BUFF(J,I,K)=ZERO
          ENDDO
        ENDDO
      ENDDO
      IF(ISSCURRENT.NE.0) GOTO 704
C
C--RUN THROUGH EVERY CELL IN THE GRID
      KT=0
      DO K=1,NLAY
        LC=LAYCON(K)
        IF(LC.EQ.3 .OR. LC.EQ.2) KT=KT+1
        DO I=1,NROW
          DO J=1,NCOL
C
C--CALCULATE FLOW FROM STORAGE (VARIABLE HEAD CELLS ONLY)
            IF(IBOUND(J,I,K).GT.0) THEN
              HSING=HNEW(J,I,K)
              IF(LC.NE.3 .AND. LC.NE.2) THEN
                RHO=SC1(J,I,K)*TLED
                STRG=RHO*HOLD(J,I,K) - RHO*HSING
              ELSE
                TP=BOTM(J,I,LBOTM(K)-1)
                RHO2=SC2(J,I,KT)*TLED
                RHO1=SC1(J,I,K)*TLED
                SOLD=RHO2
                IF(HOLD(J,I,K).GT.TP) SOLD=RHO1
                SNEW=RHO2
                IF(HSING.GT.TP) SNEW=RHO1
                STRG=SOLD*(HOLD(J,I,K)-TP) + SNEW*TP - SNEW*HSING
              ENDIF
              BUFF(J,I,K)=STRG
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--RECORD CONTENTS OF BUFFER.
  704 IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
        WRITE(IUMT3D) BUFF
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
        WRITE(IUMT3D,*) BUFF
      ENDIF
C
  705 CONTINUE
C
C--CALCULATE FLOW INTO OR OUT OF CONSTANT-HEAD CELLS
      TEXT='CNH'
      NCNH=0
C
C--CLEAR BUFFER
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            BUFF(J,I,K)=ZERO
          ENDDO
        ENDDO
      ENDDO
C
C--FOR EACH CELL IF IT IS CONSTANT HEAD COMPUTE FLOW ACROSS 6
C--FACES.
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
C
C--IF CELL IS NOT CONSTANT HEAD SKIP IT & GO ON TO NEXT CELL.
            IF(IBOUND(J,I,K).GE.0) CYCLE
            NCNH=NCNH+1
C
C--CLEAR FIELDS FOR SIX FLOW RATES.
            X1=ZERO
            X2=ZERO
            X3=ZERO
            X4=ZERO
            X5=ZERO
            X6=ZERO
C
C--CALCULATE FLOW THROUGH THE LEFT FACE
C
C--IF THERE IS AN INACTIVE CELL ON THE OTHER SIDE OF THIS
C--FACE THEN GO ON TO THE NEXT FACE.
            IF(J.EQ.1) GO TO 30
            IF(IBOUND(J-1,I,K).EQ.0) GO TO 30
            HDIFF=HNEW(J,I,K)-HNEW(J-1,I,K)
C
C--CALCULATE FLOW THROUGH THIS FACE INTO THE ADJACENT CELL.
            X1=HDIFF*CR(J-1,I,K)
C
C--CALCULATE FLOW THROUGH THE RIGHT FACE
   30       IF(J.EQ.NCOL) GO TO 60
            IF(IBOUND(J+1,I,K).EQ.0) GO TO 60
            HDIFF=HNEW(J,I,K)-HNEW(J+1,I,K)
            X2=HDIFF*CR(J,I,K)
C
C--CALCULATE FLOW THROUGH THE BACK FACE.
   60       IF(I.EQ.1) GO TO 90
            IF (IBOUND(J,I-1,K).EQ.0) GO TO 90
            HDIFF=HNEW(J,I,K)-HNEW(J,I-1,K)
            X3=HDIFF*CC(J,I-1,K)
C
C--CALCULATE FLOW THROUGH THE FRONT FACE.
   90       IF(I.EQ.NROW) GO TO 120
            IF(IBOUND(J,I+1,K).EQ.0) GO TO 120
            HDIFF=HNEW(J,I,K)-HNEW(J,I+1,K)
            X4=HDIFF*CC(J,I,K)
C
C--CALCULATE FLOW THROUGH THE UPPER FACE
  120       IF(K.EQ.1) GO TO 150
            IF (IBOUND(J,I,K-1).EQ.0) GO TO 150
            HD=HNEW(J,I,K)
            IF(LAYCON(K).NE.3 .AND. LAYCON(K).NE.2) GO TO 122
            TMP=HD
            IF(TMP.LT.BOTM(J,I,LBOTM(K)-1))
     &       HD=BOTM(J,I,LBOTM(K)-1)
  122       HDIFF=HD-HNEW(J,I,K-1)
            X5=HDIFF*CV(J,I,K-1)
C
C--CALCULATE FLOW THROUGH THE LOWER FACE.
  150       IF(K.EQ.NLAY) GO TO 180
            IF(IBOUND(J,I,K+1).EQ.0) GO TO 180
            HD=HNEW(J,I,K+1)
            IF(LAYCON(K+1).NE.3 .AND. LAYCON(K+1).NE.2) GO TO 152
            TMP=HD
            IF(TMP.LT.BOTM(J,I,LBOTM(K+1)-1))
     &       HD=BOTM(J,I,LBOTM(K+1)-1)
  152       HDIFF=HNEW(J,I,K)-HD
            X6=HDIFF*CV(J,I,K)
C
C--SUM UP FLOWS THROUGH SIX SIDES OF CONSTANT HEAD CELL.
  180       BUFF(J,I,K)=X1+X2+X3+X4+X5+X6
          ENDDO
        ENDDO
      ENDDO
C
C--RECORD CONTENTS OF BUFFER.
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT,NCNH
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT,NCNH
      ENDIF
C
C--IF THERE ARE NO CONSTANT-HEAD CELLS THEN SKIP
      IF(NCNH.LE.0) GOTO 1000
C
C--WRITE CONSTANT-HEAD CELL LOCATIONS AND RATES
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            IF(IBOUND(J,I,K).GE.0) CYCLE
            IF(ILMTFMT.EQ.0) THEN
              WRITE(IUMT3D)   K,I,J,BUFF(J,I,K)            
            ELSEIF(ILMTFMT.EQ.1) THEN
              WRITE(IUMT3D,*) K,I,J,BUFF(J,I,K)              
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--RETURN
 1000 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE LMT7LPF7(ILMTFMT,ISSMT3D,IUMT3D,KSTP,KPER,IGRID)
C *********************************************************************
C SAVE FLOW ACROSS THREE CELL INTERFACES (QXX, QYY, QZZ), FLOW RATE TO
C OR FROM TRANSIENT FLUID-STORAGE (QSTO), AND LOCATIONS AND FLOW RATES
C OF CONSTANT-HEAD CELLS FOR USE BY MT3D.  THIS SUBROUTINE IS CALLED
C ONLY IF THE 'LPF' PACKAGE IS USED IN MODFLOW.
C *********************************************************************
C Modified from Harbaugh(2005)
C last modified: 08-08-2008
C
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,ISSFLG,IBOUND,HNEW,HOLD,
     &                      BUFF,CR,CC,CV,BOTM,LBOTM
      USE GWFBASMODULE,ONLY:DELT
      USE GWFLPFMODULE,ONLY:LAYTYP,SC1,SC2
      CHARACTER*16 TEXT
      DOUBLE PRECISION HD
C
C--SET POINTERS FOR THE CURRENT GRID      
cswm: already set GWF2LPF7BDS      CALL SGWF2LPF7PNT(IGRID)
C      
C--GET STEADY-STATE FLAG FOR THE CURRENT STRESS PERIOD
      ISSCURRENT=ISSFLG(KPER)      
C
C--CALCULATE AND SAVE SATURATED THICKNESS
      TEXT='THKSAT'
      ZERO=0.
      ONE=1.
C
C--INITIALIZE BUFF ARRAY WITH 1.E30 FOR INACTIVE CELLS
C--OR FLAG -111 FOR ACTIVE CELLS
      FlagInactive=1.E30
      FlagActive=-111.
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            IF(IBOUND(J,I,K).EQ.0) THEN
              BUFF(J,I,K)=FlagInactive
            ELSE
              BUFF(J,I,K)=FlagActive
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--CALCULATE SATURATED THICKNESS FOR UNCONFINED/CONVERTIBLE
C--LAYERS AND STORE IN ARRAY BUFF
      DO K=1,NLAY
        IF(LAYTYP(K).EQ.0) CYCLE
        DO I=1,NROW
          DO J=1,NCOL
            IF(IBOUND(J,I,K).NE.0) THEN
              TMP=HNEW(J,I,K)
              BUFF(J,I,K)=TMP-BOTM(J,I,LBOTM(K))
              THKLAY=BOTM(J,I,LBOTM(K)-1)-BOTM(J,I,LBOTM(K))
              IF(BUFF(J,I,K).GT.THKLAY) BUFF(J,I,K)=THKLAY
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--SAVE THE CONTENTS OF THE BUFFER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
        WRITE(IUMT3D) BUFF
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
        WRITE(IUMT3D,*) BUFF
      ENDIF
C
C--CALCULATE AND SAVE FLOW ACROSS RIGHT FACE
      NCM1=NCOL-1
      IF(NCM1.LT.1) GO TO 405
      TEXT='QXX'
C
C--CLEAR THE BUFFER
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            BUFF(J,I,K)=ZERO
          ENDDO
        ENDDO
      ENDDO
C
C--FOR EACH CELL
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCM1
            IF(IBOUND(J,I,K).NE.0.AND.IBOUND(J+1,I,K).NE.0) THEN
              HDIFF=HNEW(J,I,K)-HNEW(J+1,I,K)
              BUFF(J,I,K)=HDIFF*CR(J,I,K)
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--RECORD CONTENTS OF BUFFER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
        WRITE(IUMT3D) BUFF
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
        WRITE(IUMT3D,*) BUFF
      ENDIF
C
  405 CONTINUE
C
C--CALCULATE AND SAVE FLOW ACROSS FRONT FACE
      NRM1=NROW-1
      IF(NRM1.LT.1) GO TO 505
      TEXT='QYY'
C
C--CLEAR THE BUFFER
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            BUFF(J,I,K)=ZERO
          ENDDO
        ENDDO
      ENDDO
C
C--FOR EACH CELL
      DO K=1,NLAY
        DO I=1,NRM1
          DO J=1,NCOL
            IF(IBOUND(J,I,K).NE.0.AND.IBOUND(J,I+1,K).NE.0) THEN
              HDIFF=HNEW(J,I,K)-HNEW(J,I+1,K)
              BUFF(J,I,K)=HDIFF*CC(J,I,K)
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--RECORD CONTENTS OF BUFFER.
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
        WRITE(IUMT3D) BUFF
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
        WRITE(IUMT3D,*) BUFF
      ENDIF
C
  505 CONTINUE
C
C--CALCULATE AND SAVE FLOW ACROSS FRONT FACE
      NLM1=NLAY-1
      IF(NLM1.LT.1) GO TO 700
      TEXT='QZZ'
C
C--CLEAR THE BUFFER
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            BUFF(J,I,K)=ZERO
          ENDDO
        ENDDO
      ENDDO
C
C--FOR EACH CELL CALCULATE FLOW THRU LOWER FACE & STORE IN BUFFER
      DO K=1,NLM1
        DO I=1,NROW
          DO J=1,NCOL
            IF(IBOUND(J,I,K).NE.0.AND.IBOUND(J,I,K+1).NE.0) THEN
              HD=HNEW(J,I,K+1)
              IF(LAYTYP(K+1).NE.0) THEN
                TMP=HD
                TOP=BOTM(J,I,LBOTM(K+1)-1)
                IF(TMP.LT.TOP) HD=TOP
              ENDIF
              HDIFF=HNEW(J,I,K)-HD
              BUFF(J,I,K)=HDIFF*CV(J,I,K)
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--RECORD CONTENTS OF BUFFER.
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
        WRITE(IUMT3D) BUFF
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
        WRITE(IUMT3D,*) BUFF
      ENDIF
C
  700 CONTINUE
C
C--CALCULATE AND SAVE GROUNDWATER STORAGE IF TRANSIENT
      IF(ISSMT3D.NE.0) GO TO 705
      TEXT='STO'
C
C--INITIALIZE AND CLEAR BUFFER           
cswm: moved below      TLED=ONE/DELT
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            BUFF(J,I,K)=ZERO
          ENDDO
        ENDDO
      ENDDO
      IF(ISSCURRENT.NE.0) GOTO 704
C
C--RUN THROUGH EVERY CELL IN THE GRID
      TLED=ONE/DELT !swm: moved after check for transient sp
      KT=0
      DO K=1,NLAY
        LC=LAYTYP(K)
        IF(LC.NE.0) KT=KT+1
        DO I=1,NROW
          DO J=1,NCOL
C
C--CALCULATE FLOW FROM STORAGE (VARIABLE HEAD CELLS ONLY)
            IF(IBOUND(J,I,K).GT.0) THEN
              HSING=HNEW(J,I,K)
              IF(LC.EQ.0) THEN
                RHO=SC1(J,I,K)*TLED
                STRG=RHO*HOLD(J,I,K) - RHO*HSING
              ELSE
                TP=BOTM(J,I,LBOTM(K)-1)
                RHO2=SC2(J,I,KT)*TLED
                RHO1=SC1(J,I,K)*TLED
                SOLD=RHO2
                IF(HOLD(J,I,K).GT.TP) SOLD=RHO1
                SNEW=RHO2
                IF(HSING.GT.TP) SNEW=RHO1
                STRG=SOLD*(HOLD(J,I,K)-TP) + SNEW*TP - SNEW*HSING
              ENDIF
              BUFF(J,I,K)=STRG
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--RECORD CONTENTS OF BUFFER.
  704 IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
        WRITE(IUMT3D) BUFF
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
        WRITE(IUMT3D,*) BUFF
      ENDIF
C
  705 CONTINUE
C
C--CALCULATE FLOW INTO OR OUT OF CONSTANT-HEAD CELLS
      TEXT='CNH'
      NCNH=0
C
C--CLEAR BUFFER
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            BUFF(J,I,K)=ZERO
          ENDDO
        ENDDO
      ENDDO
C
C--FOR EACH CELL IF IT IS CONSTANT HEAD COMPUTE FLOW ACROSS 6
C--FACES.
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
C
C--IF CELL IS NOT CONSTANT HEAD SKIP IT & GO ON TO NEXT CELL.
            IF(IBOUND(J,I,K).GE.0) CYCLE
            NCNH=NCNH+1
C
C--CLEAR FIELDS FOR SIX FLOW RATES.
            X1=ZERO
            X2=ZERO
            X3=ZERO
            X4=ZERO
            X5=ZERO
            X6=ZERO
C
C--CALCULATE FLOW THROUGH THE LEFT FACE
C
C--IF THERE IS AN INACTIVE CELL ON THE OTHER SIDE OF THIS
C--FACE THEN GO ON TO THE NEXT FACE.
            IF(J.EQ.1) GO TO 30
            IF(IBOUND(J-1,I,K).EQ.0) GO TO 30
            HDIFF=HNEW(J,I,K)-HNEW(J-1,I,K)
C
C--CALCULATE FLOW THROUGH THIS FACE INTO THE ADJACENT CELL.
            X1=HDIFF*CR(J-1,I,K)
C
C--CALCULATE FLOW THROUGH THE RIGHT FACE
   30       IF(J.EQ.NCOL) GO TO 60
            IF(IBOUND(J+1,I,K).EQ.0) GO TO 60
            HDIFF=HNEW(J,I,K)-HNEW(J+1,I,K)
            X2=HDIFF*CR(J,I,K)
C
C--CALCULATE FLOW THROUGH THE BACK FACE.
   60       IF(I.EQ.1) GO TO 90
            IF (IBOUND(J,I-1,K).EQ.0) GO TO 90
            HDIFF=HNEW(J,I,K)-HNEW(J,I-1,K)
            X3=HDIFF*CC(J,I-1,K)
C
C--CALCULATE FLOW THROUGH THE FRONT FACE.
   90       IF(I.EQ.NROW) GO TO 120
            IF(IBOUND(J,I+1,K).EQ.0) GO TO 120
            HDIFF=HNEW(J,I,K)-HNEW(J,I+1,K)
            X4=HDIFF*CC(J,I,K)
C
C--CALCULATE FLOW THROUGH THE UPPER FACE
  120       IF(K.EQ.1) GO TO 150
            IF (IBOUND(J,I,K-1).EQ.0) GO TO 150
            HD=HNEW(J,I,K)
            IF(LAYTYP(K).EQ.0) GO TO 122
            TMP=HD
            TOP=BOTM(J,I,LBOTM(K)-1)
            IF(TMP.LT.TOP) HD=TOP
  122       HDIFF=HD-HNEW(J,I,K-1)
            X5=HDIFF*CV(J,I,K-1)
C
C--CALCULATE FLOW THROUGH THE LOWER FACE.
  150       IF(K.EQ.NLAY) GO TO 180
            IF(IBOUND(J,I,K+1).EQ.0) GO TO 180
            HD=HNEW(J,I,K+1)
            IF(LAYTYP(K+1).EQ.0) GO TO 152
            TMP=HD
            TOP=BOTM(J,I,LBOTM(K+1)-1)
            IF(TMP.LT.TOP) HD=TOP
  152       HDIFF=HNEW(J,I,K)-HD
            X6=HDIFF*CV(J,I,K)
C
C--SUM UP FLOWS THROUGH SIX SIDES OF CONSTANT HEAD CELL.
  180       BUFF(J,I,K)=X1+X2+X3+X4+X5+X6
          ENDDO
        ENDDO
      ENDDO
C
C--RECORD CONTENTS OF BUFFER.
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT,NCNH
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT,NCNH
      ENDIF
C
C--IF THERE ARE NO CONSTANT-HEAD CELLS THEN SKIP
      IF(NCNH.LE.0) GOTO 1000
C
C--WRITE CONSTANT-HEAD CELL LOCATIONS AND RATES
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            IF(IBOUND(J,I,K).GE.0) CYCLE
            IF(ILMTFMT.EQ.0) THEN
              WRITE(IUMT3D)   K,I,J,BUFF(J,I,K)
            ELSEIF(ILMTFMT.EQ.1) THEN
              WRITE(IUMT3D,*) K,I,J,BUFF(J,I,K)
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--RETURN
 1000 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE LMT7HUF7(ILMTFMT,ISSMT3D,IUMT3D,KSTP,KPER,ILVDA,IGRID)
C **********************************************************************
C SAVE FLOW ACROSS THREE CELL INTERFACES (QXX, QYY, QZZ), FLOW RATE TO
C OR FROM TRANSIENT FLUID-STORAGE (QSTO), AND LOCATIONS AND FLOW RATES
C OF CONSTANT-HEAD CELLS FOR USE BY MT3D.  THIS SUBROUTINE IS CALLED
C ONLY IF THE 'HUF' PACKAGE IS USED IN MODFLOW.
C **********************************************************************
C Modified from Anderman and Hill (2000), Harbaugh (2005)
C last modified: 08-08-2008
C
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,ISSFLG,IBOUND,HNEW,HOLD,BOTM,
     &                      LBOTM,DELR,DELC,BUFF,IOUT,CR,CC,CV
      USE GWFBASMODULE,ONLY:DELT
      USE GWFHUFMODULE,ONLY:LTHUF,SC1,HUFTHK,NHUF,VDHT      
      CHARACTER*16 TEXT
      DOUBLE PRECISION HD,DFL,DFR,DFT,DFB,HN
C    
C--SET POINTERS FOR THE CURRENT GRID   
cswm: already set in GWF2HUF7BDS      CALL SGWF2HUF7PNT(IGRID) 
C      
C--GET STEADY-STATE FLAG FOR THE CURRENT STRESS PERIOD      
      ISSCURRENT=ISSFLG(KPER)
C
C--CALCULATE AND SAVE SATURATED THICKNESS
      TEXT='THKSAT'
      ZERO=0.
      ONE=1.   
C
C--INITIALIZE BUFF ARRAY WITH 1.E30 FOR INACTIVE CELLS
C--OR FLAG -111 FOR ACTIVE CELLS
      FlagInactive=1.E30
      FlagActive=-111.   
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            IF(IBOUND(J,I,K).EQ.0) THEN
              BUFF(J,I,K)=FlagInactive
            ELSE
              BUFF(J,I,K)=FlagActive
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--CALCULATE SATURATED THICKNESS FOR UNCONFINED/CONVERTIBLE
C--LAYERS AND STORE IN ARRAY BUFF
      DO K=1,NLAY
        IF(LTHUF(K).EQ.0) CYCLE
        DO I=1,NROW
          DO J=1,NCOL
            IF(IBOUND(J,I,K).NE.0) THEN
              TMP=HNEW(J,I,K)
              BUFF(J,I,K)=TMP-BOTM(J,I,LBOTM(K))
              THKLAY=BOTM(J,I,LBOTM(K)-1)-BOTM(J,I,LBOTM(K))
              IF(BUFF(J,I,K).GT.THKLAY) BUFF(J,I,K)=THKLAY
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--SAVE THE CONTENTS OF THE BUFFER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
        WRITE(IUMT3D) BUFF
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
        WRITE(IUMT3D,*) BUFF
      ENDIF
C
C--CALCULATE AND SAVE FLOW ACROSS RIGHT FACE
      NCM1=NCOL-1
      IF(NCM1.LT.1) GO TO 405
      TEXT='QXX'
C
C--CLEAR THE BUFFER
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            BUFF(J,I,K)=ZERO
          ENDDO
        ENDDO
      ENDDO
C
C--FOR EACH CELL
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCM1
            IF(IBOUND(J,I,K).NE.0.AND.IBOUND(J+1,I,K).NE.0) THEN        
              if(ILVDA.gt.0) then
                CALL SGWF2HUF7VDF9(I,J,K,VDHT,HNEW,IBOUND,
     &           NLAY,NROW,NCOL,DFL,DFR,DFT,DFB)
                BUFF(J,I,K) = DFR
              else                       
                HDIFF=HNEW(J,I,K)-HNEW(J+1,I,K)
                BUFF(J,I,K)=HDIFF*CR(J,I,K)
              endif  
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--RECORD CONTENTS OF BUFFER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
        WRITE(IUMT3D) BUFF
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
        WRITE(IUMT3D,*) BUFF
      ENDIF
C
  405 CONTINUE
C
C--CALCULATE AND SAVE FLOW ACROSS FRONT FACE
      NRM1=NROW-1
      IF(NRM1.LT.1) GO TO 505
      TEXT='QYY'
C
C--CLEAR THE BUFFER
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            BUFF(J,I,K)=ZERO
          ENDDO
        ENDDO
      ENDDO
C
C--FOR EACH CELL
      DO K=1,NLAY
        DO I=1,NRM1
          DO J=1,NCOL
            IF(IBOUND(J,I,K).NE.0.AND.IBOUND(J,I+1,K).NE.0) THEN
              if(ILVDA.gt.0) then
                CALL SGWF2HUF7VDF9(I,J,K,VDHT,HNEW,IBOUND,
     &           NLAY,NROW,NCOL,DFL,DFR,DFT,DFB)
                BUFF(J,I,K) = DFT
              else                        
                HDIFF=HNEW(J,I,K)-HNEW(J,I+1,K)
                BUFF(J,I,K)=HDIFF*CC(J,I,K)
              endif  
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--RECORD CONTENTS OF BUFFER.
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
        WRITE(IUMT3D) BUFF
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
        WRITE(IUMT3D,*) BUFF
      ENDIF
C
  505 CONTINUE
C
C--CALCULATE AND SAVE FLOW ACROSS LOWER FACE
      NLM1=NLAY-1
      IF(NLM1.LT.1) GO TO 700
      TEXT='QZZ'
C
C--CLEAR THE BUFFER
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            BUFF(J,I,K)=ZERO
          ENDDO
        ENDDO
      ENDDO
C
C--FOR EACH CELL 
      DO K=1,NLM1
        DO I=1,NROW
          DO J=1,NCOL
            IF(IBOUND(J,I,K).NE.0.AND.IBOUND(J,I,K+1).NE.0) THEN
              HD=HNEW(J,I,K+1)
              IF(LTHUF(K+1).NE.0) THEN
                TMP=HD
                TOP=BOTM(J,I,LBOTM(K+1)-1)
                IF(TMP.LT.TOP) HD=TOP
              ENDIF
              HDIFF=HNEW(J,I,K)-HD
              BUFF(J,I,K)=HDIFF*CV(J,I,K)
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--RECORD CONTENTS OF BUFFER.
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
        WRITE(IUMT3D) BUFF
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
        WRITE(IUMT3D,*) BUFF
      ENDIF
C
  700 CONTINUE
C
C--CALCULATE AND SAVE GROUNDWATER STORAGE IF TRANSIENT
      IF(ISSMT3D.NE.0) GO TO 705
      TEXT='STO'
C
C--INITIALIZE and CLEAR BUFFER
      TLED=ONE/DELT
      DO K=1,NLAY 
        DO I=1,NROW
          DO J=1,NCOL
            BUFF(J,I,K)=ZERO
          ENDDO
        ENDDO
      ENDDO
      IF(ISSCURRENT.NE.0) GOTO 704
C
C5------LOOP THROUGH EVERY CELL IN THE GRID.
      KT=0
      DO K=1,NLAY
        LC=LTHUF(K)
        IF(LC.NE.0) KT=KT+1
        DO I=1,NROW
          DO J=1,NCOL
C
C6------SKIP NO-FLOW AND CONSTANT-HEAD CELLS.
            IF(IBOUND(J,I,K).LE.0) CYCLE
            HN=HNEW(J,I,K)
            HO=HOLD(J,I,K)
            STRG=ZERO
C
C7-----CHECK LAYER TYPE TO SEE IF ONE STORAGE CAPACITY OR TWO.
            IF(LC.EQ.0) GO TO 285
            TOP=BOTM(J,I,LBOTM(K)-1)
            BOT=BOTM(J,I,LBOTM(K))
            IF(HO.GT.TOP.AND.HN.GT.TOP) GOTO 285
C
C7A----TWO STORAGE CAPACITIES.
C---------------Compute SC1 Component
            IF(HO.GT.TOP) THEN
              STRG=SC1(J,I,K)*(HO-TOP)*TLED
            ELSEIF(HN.GT.TOP) THEN
              STRG=SC1(J,I,K)*TLED*(TOP-HN)
            ENDIF
C---------------Compute SC2 Component
            CALL SGWF2HUF7SC2(1,J,I,K,TOP,BOT,HN,HO,TLED,CHCOF,STRG,
     &       HUFTHK,NCOL,NROW,NHUF,DELR(J)*DELC(I),IOUT)          
C------STRG=SOLD*(HOLD(J,I,K)-TP) + SNEW*TP - SNEW*HSING
            GOTO 288
C
C7B----ONE STORAGE CAPACITY.
  285       RHO=SC1(J,I,K)*TLED
            STRG=RHO*(HO-HN)
C
C8-----STORE CELL-BY-CELL FLOW IN BUFFER
  288       BUFF(J,I,K)=STRG
C
          ENDDO
        ENDDO
      ENDDO
C
C--RECORD CONTENTS OF BUFFER.
  704 IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
        WRITE(IUMT3D) BUFF
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
        WRITE(IUMT3D,*) BUFF
      ENDIF
C
  705 CONTINUE
C
C--CALCULATE FLOW INTO OR OUT OF CONSTANT-HEAD CELLS
      TEXT='CNH'
      NCNH=0
C
C--CLEAR BUFFER
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            BUFF(J,I,K)=ZERO
          ENDDO
        ENDDO
      ENDDO
C
C--FOR EACH CELL IF IT IS CONSTANT HEAD COMPUTE FLOW ACROSS 6
C--FACES.
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
C
C--IF CELL IS NOT CONSTANT HEAD SKIP IT & GO ON TO NEXT CELL.
            IF (IBOUND(J,I,K).GE.0) CYCLE
            NCNH=NCNH+1
C
C--CLEAR FIELDS FOR SIX FLOW RATES.
            X1=ZERO
            X2=ZERO
            X3=ZERO
            X4=ZERO
            X5=ZERO
            X6=ZERO
C            
C--COMPUTE HORIZONTAL FLUXES IF THE LVDA CAPABILITY IS USED            
            if(ILVDA.gt.0)
     &       CALL SGWF2HUF7VDF9(I,J,K,VDHT,HNEW,IBOUND,
     &       NLAY,NROW,NCOL,DFL,DFR,DFT,DFB)                        
C
C--CALCULATE FLOW THROUGH THE LEFT FACE
C
C--IF THERE IS AN INACTIVE CELL ON THE OTHER SIDE OF THIS
C--FACE THEN GO ON TO THE NEXT FACE.
            IF(J.EQ.1) GO TO 30
            IF(IBOUND(J-1,I,K).EQ.0) GO TO 30
C
C--CALCULATE FLOW THROUGH THIS FACE INTO THE ADJACENT CELL.
            if(ILVDA.gt.0) then
              X1 = -DFL
            else
              HDIFF=HNEW(J,I,K)-HNEW(J-1,I,K)            
              X1=HDIFF*CR(J-1,I,K)
            endif  
C
C--CALCULATE FLOW THROUGH THE RIGHT FACE
   30       IF(J.EQ.NCOL) GO TO 60
            IF(IBOUND(J+1,I,K).EQ.0) GO TO 60
            if(ILVDA.gt.0) then
              X2 = DFR
            else                       
              HDIFF=HNEW(J,I,K)-HNEW(J+1,I,K)
              X2=HDIFF*CR(J,I,K)
            endif  
C
C--CALCULATE FLOW THROUGH THE BACK FACE.
   60       IF(I.EQ.1) GO TO 90
            IF (IBOUND(J,I-1,K).EQ.0) GO TO 90
            if(ILVDA.gt.0) then
              X3 = -DFT
            else                       
              HDIFF=HNEW(J,I,K)-HNEW(J,I-1,K)
              X3=HDIFF*CC(J,I-1,K)
            endif  
C
C--CALCULATE FLOW THROUGH THE FRONT FACE.
   90       IF(I.EQ.NROW) GO TO 120
            IF(IBOUND(J,I+1,K).EQ.0) GO TO 120
            if(ILVDA.gt.0) then
              X4 = DFB
            else             
              HDIFF=HNEW(J,I,K)-HNEW(J,I+1,K)
              X4=HDIFF*CC(J,I,K)
            endif  
C
C--CALCULATE FLOW THROUGH THE UPPER FACE
  120       IF(K.EQ.1) GO TO 150
            IF (IBOUND(J,I,K-1).EQ.0) GO TO 150
            HD=HNEW(J,I,K)
            IF(LTHUF(K).EQ.0) GO TO 122
            TMP=HD
            TOP=BOTM(J,I,LBOTM(K)-1)
            IF(TMP.LT.TOP) HD=TOP
  122       HDIFF=HD-HNEW(J,I,K-1)
            X5=HDIFF*CV(J,I,K-1)
C
C--CALCULATE FLOW THROUGH THE LOWER FACE.
  150       IF(K.EQ.NLAY) GO TO 180
            IF(IBOUND(J,I,K+1).EQ.0) GO TO 180
            HD=HNEW(J,I,K+1)
            IF(LTHUF(K+1).EQ.0) GO TO 152
            TMP=HD
            TOP=BOTM(J,I,LBOTM(K+1)-1)
            IF(TMP.LT.TOP) HD=TOP
  152       HDIFF=HNEW(J,I,K)-HD
            X6=HDIFF*CV(J,I,K)
C
C--SUM UP FLOWS THROUGH SIX SIDES OF CONSTANT HEAD CELL.
  180       BUFF(J,I,K)=X1+X2+X3+X4+X5+X6
C
          ENDDO
        ENDDO
      ENDDO
C
C--RECORD CONTENTS OF BUFFER.
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT,NCNH
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT,NCNH
      ENDIF
C
C--IF THERE ARE NO CONSTANT-HEAD CELLS THEN SKIP
      IF(NCNH.LE.0) GOTO 1000
C
C--WRITE CONSTANT-HEAD CELL LOCATIONS AND RATES
      DO K=1,NLAY
        DO I=1,NROW
          DO J=1,NCOL
            IF(IBOUND(J,I,K).GE.0) CYCLE
            IF(ILMTFMT.EQ.0) THEN
              WRITE(IUMT3D)   K,I,J,BUFF(J,I,K)
            ELSEIF(ILMTFMT.EQ.1) THEN
              WRITE(IUMT3D,*) K,I,J,BUFF(J,I,K)
            ENDIF              
          ENDDO
        ENDDO
      ENDDO
C
C--RETURN
 1000 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE LMT7WEL7(ILMTFMT,IUMT3D,KSTP,KPER,IGRID)
C *********************************************************************
C SAVE WELL CELL LOCATIONS AND VOLUMETRIC FLOW RATES FOR USE BY MT3D.
C *********************************************************************
C Modified from  Harbaugh (2005)
C last modified: 08-08-2008
C
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,IBOUND
      USE GWFWELMODULE,ONLY:NWELLS,WELL      
      CHARACTER*16 TEXT
C      
C--SET POINTERS FOR THE CURRENT GRID   
cswm: already set in      CALL SGWF2WEL7PNT(IGRID)
C      
      TEXT='WEL'   
      ZERO=0.
C
C--WRITE AN IDENTIFYING HEADER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT,NWELLS
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT,NWELLS
      ENDIF
C
C--IF THERE ARE NO WELLS RETURN
      IF(NWELLS.LE.0) GO TO 9999
C
C--WRITE WELL LOCATION AND RATE ONE AT A TIME
      DO L=1,NWELLS
        IL=WELL(1,L)
        IR=WELL(2,L)
        IC=WELL(3,L)
C
C--IF CELL IS EXTERNAL Q=0
        Q=ZERO
        IF(IBOUND(IC,IR,IL).GT.0) Q=WELL(4,L)
        IF(ILMTFMT.EQ.0) THEN
          WRITE(IUMT3D) IL,IR,IC,Q
        ELSEIF(ILMTFMT.EQ.1) THEN
          WRITE(IUMT3D,*) IL,IR,IC,Q
        ENDIF
      ENDDO
C
C--RETURN
 9999 RETURN
      END
C
C
      SUBROUTINE LMT7DRN7(ILMTFMT,IUMT3D,KSTP,KPER,IGRID)
C ********************************************************************
C SAVE DRAIN CELL LOCATIONS AND VOLUMETRIC FLOW RATES FOR USE BY MT3D.
C ********************************************************************
C Modified from Harbaugh (2005)
C last modified: 08-08-2008
C
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,IBOUND,HNEW,BUFF
      USE GWFDRNMODULE,ONLY:NDRAIN,DRAI
      CHARACTER*16 TEXT
      DOUBLE PRECISION HHNEW,EEL,CCDRN,CEL,QQ
C    
C--SET POINTERS FOR THE CURRENT GRID
c swm: already set in GWF2DRN7BD      CALL SGWF2DRN7PNT(IGRID)
C      
      TEXT='DRN'
      ZERO=0.
C
C--WRITE AN IDENTIFYING HEADER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT,NDRAIN
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT,NDRAIN
      ENDIF
C
C--IF THERE ARE NO DRAINS THEN SKIP
      IF(NDRAIN.LE.0) GO TO 9999
C
C--FOR EACH DRAIN ACCUMULATE DRAIN FLOW
      DO L=1,NDRAIN
C
C--GET LAYER, ROW & COLUMN OF CELL CONTAINING REACH.
        IL=DRAI(1,L)
        IR=DRAI(2,L)
        IC=DRAI(3,L)
        QQ=ZERO
C
C--CALCULATE Q FOR ACTIVE CELLS
        IF(IBOUND(IC,IR,IL).GT.0) THEN
C
C--GET DRAIN PARAMETERS FROM DRAIN LIST.
          EL=DRAI(4,L)
          EEL=EL
          C=DRAI(5,L)
          CCDRN=C
          HHNEW=HNEW(IC,IR,IL)
          CEL=C*EL
C
C--IF HEAD LOWER THAN DRAIN THEN FORGET THIS CELL.
C--OTHERWISE, CALCULATE Q=C*(EL-HHNEW).
          IF(HHNEW.GT.EEL) QQ=CEL - CCDRN*HHNEW
        ENDIF
        Q=QQ
C
C--WRITE DRAIN LOCATION AND RATE
        IF(ILMTFMT.EQ.0) THEN
          WRITE(IUMT3D)   IL,IR,IC,Q
        ELSEIF(ILMTFMT.EQ.1) THEN
          WRITE(IUMT3D,*) IL,IR,IC,Q
        ENDIF  
C        
      ENDDO
C
C--RETURN
 9999 RETURN
      END
C
C
      SUBROUTINE LMT7RIV7(ILMTFMT,IUMT3D,KSTP,KPER,IGRID)
C *********************************************************************
C SAVE RIVER CELL LOCATIONS AND VOLUMETRIC FLOW RATES FOR USE BY MT3D.
C *********************************************************************
C Modified from Harbaugh (2005)
C last modified: 08-08-2008
C
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,IBOUND,HNEW
      USE GWFRIVMODULE,ONLY:NRIVER,RIVR
      CHARACTER*16 TEXT      
      DOUBLE PRECISION HHNEW,CHRIV,RRBOT,CCRIV
C
C--SET POINTERS FOR THE CURRENT GRID
c swm: already set in GWF2RIV7BD      CALL SGWF2RIV7PNT(IGRID)      
C      
      TEXT='RIV'      
      ZERO=0.
C
C--WRITE AN IDENTIFYING HEADER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT,NRIVER
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT,NRIVER
      ENDIF
C
C--IF NO REACHES SKIP
      IF(NRIVER.LE.0) GO TO 9999
C
C--FOR EACH RIVER REACH ACCUMULATE RIVER FLOW
      DO L=1,NRIVER
C
C--GET LAYER, ROW & COLUMN OF CELL CONTAINING REACH.
        IL=RIVR(1,L)
        IR=RIVR(2,L)
        IC=RIVR(3,L)
C
C--IF CELL IS EXTERNAL RATE=0
        IF(IBOUND(IC,IR,IL).LE.0) THEN
          RATE=ZERO
C
C--GET RIVER PARAMETERS FROM RIVER LIST.
        ELSE
          HRIV=RIVR(4,L)
          CRIV=RIVR(5,L)
          RBOT=RIVR(6,L)
          HHNEW=HNEW(IC,IR,IL)
          CHRIV=CRIV*HRIV
          CCRIV=CRIV
          RRBOT=RBOT
C
C--COMPARE HEAD IN AQUIFER TO BOTTOM OF RIVERBED.
C
C--AQUIFER HEAD > BOTTOM THEN RATE=CRIV*(HRIV-HNEW).
          IF(HHNEW.GT.RRBOT) RATE=CHRIV-CCRIV*HHNEW
C
C--AQUIFER HEAD < BOTTOM THEN RATE=CRIV*(HRIV-RBOT)
          IF(HHNEW.LE.RRBOT) RATE=CRIV*(HRIV-RBOT)
        ENDIF
C
C--WRITE RIVER REACH LOCATION AND RATE
        IF(ILMTFMT.EQ.0) THEN
          WRITE(IUMT3D) IL,IR,IC,RATE
        ELSEIF(ILMTFMT.EQ.1) THEN
          WRITE(IUMT3D,*) IL,IR,IC,RATE
        ENDIF
C        
      ENDDO
C
C--RETURN
 9999 RETURN
      END
C
C
      SUBROUTINE LMT7RCH7(ILMTFMT,IUMT3D,KSTP,KPER,IGRID)
C *******************************************************************
C SAVE RECHARGE LAYER LOCATION AND VOLUMETRIC FLOW RATES
C FOR USE BY MT3D.
C *******************************************************************
C Modified from Harbaugh (2005)
C last modified: 08-08-2008
C
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,IBOUND,BUFF
      USE GWFRCHMODULE,ONLY:NRCHOP,RECH,IRCH
      CHARACTER*16 TEXT
C
C--SET POINTERS FOR THE CURRENT GRID
c swm: already set in GWF2RCH7BD      CALL SGWF2RCH7PNT(IGRID)   
C         
      TEXT='RCH'
      ZERO=0.
C
C--WRITE AN IDENTIFYING HEADER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
      ENDIF
C
C--CLEAR THE BUFFER.
      DO IL=1,NLAY
        DO IR=1,NROW
          DO IC=1,NCOL
            BUFF(IC,IR,IL)=ZERO
          ENDDO
        ENDDO
      ENDDO
C
C--IF NRCHOP=1 RECH GOES INTO LAYER 1.
      IF(NRCHOP.EQ.1) THEN
        IL=1
        IF(ILMTFMT.EQ.0) WRITE(IUMT3D)   ((IL,J=1,NCOL),I=1,NROW)
        IF(ILMTFMT.EQ.1) WRITE(IUMT3D,*) ((IL,J=1,NCOL),I=1,NROW)
C
C--STORE RECH RATE IN BUFF FOR ACTIVE CELLS
        DO I=1,NROW
          DO J=1,NCOL
            IF(IBOUND(J,I,1).GT.0) BUFF(J,I,1)=RECH(J,I)
          ENDDO
        ENDDO
        IF(ILMTFMT.EQ.0) THEN 
          WRITE(IUMT3D)   ((BUFF(J,I,1),J=1,NCOL),I=1,NROW)
        ELSEIF(ILMTFMT.EQ.1) THEN
          WRITE(IUMT3D,*) ((BUFF(J,I,1),J=1,NCOL),I=1,NROW)
        ENDIF
C
C--IF NRCHOP=2 OR 3 RECH IS IN LAYER SHOWN IN INDICATOR ARRAY(IRCH).
      ELSEIF(NRCHOP.NE.1) THEN
        IF(ILMTFMT.EQ.0) THEN
          WRITE(IUMT3D)   ((IRCH(J,I),J=1,NCOL),I=1,NROW)
        ELSEIF(ILMTFMT.EQ.1) THEN
          WRITE(IUMT3D,*) ((IRCH(J,I),J=1,NCOL),I=1,NROW)
        ENDIF  
C
C--STORE RECH RATE IN BUFF FOR ACTIVE CELLS
        DO I=1,NROW
          DO J=1,NCOL
            IL=IRCH(J,I)
            IF(IL.EQ.0) CYCLE
            IF(IBOUND(J,I,IL).GT.0) THEN
              BUFF(J,I,1)=RECH(J,I)
            ENDIF
          ENDDO
        ENDDO
        IF(ILMTFMT.EQ.0) THEN 
          WRITE(IUMT3D)   ((BUFF(J,I,1),J=1,NCOL),I=1,NROW)
        ELSEIF(ILMTFMT.EQ.1) THEN
          WRITE(IUMT3D,*) ((BUFF(J,I,1),J=1,NCOL),I=1,NROW)
        ENDIF  
      ENDIF
C
C--RETURN
      RETURN
      END
C
C
      SUBROUTINE LMT7EVT7(ILMTFMT,IUMT3D,KSTP,KPER,IGRID)
C ******************************************************************
C SAVE EVAPOTRANSPIRATION LAYER LOCATION AND VOLUMETRIC FLOW RATES
C FOR USE BY MT3D.
C ******************************************************************
C Modified from Harbaugh (2005)
C last modified: 08-08-2008
C
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,IBOUND,HNEW,BUFF      
      USE GWFEVTMODULE,ONLY:NEVTOP,EVTR,EXDP,SURF,IEVT
      CHARACTER*16 TEXT
      DOUBLE PRECISION QQ,HH,XX,DD,SS,HHCOF,RRHS      
C   
C--SET POINTERS FOR THE CURRENT GRID
c swm: already set in GWF2EVT7BD      CALL SGWF2EVT7PNT(IGRID)
C            
      TEXT='EVT'
      ZERO=0.      
C
C--WRITE AN IDENTIFYING HEADER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
      ENDIF
C
C--CLEAR THE BUFFER.
      DO IL=1,NLAY
        DO IR=1,NROW
          DO IC=1,NCOL
            BUFF(IC,IR,IL)=ZERO
          ENDDO
        ENDDO
      ENDDO   
C
C--PROCESS EACH HORIZONTAL CELL LOCATION
C--AND STORE ET RATES IN BUFFER (IC,IR,1)
      DO IR=1,NROW
        DO IC=1,NCOL
C
C--IF OPTION 1 SET THE LAYER INDEX EQUAL TO 1
          IF(NEVTOP.EQ.1) THEN
            IL=1
C
C--IF OPTION 2 OR 3 GET LAYER INDEX FROM IEVT ARRAY
          ELSEIF(NEVTOP.NE.1) THEN
            IL=IEVT(IC,IR)
            IF(IL.EQ.0) CYCLE
          ENDIF
C
C--IF CELL IS EXTERNAL THEN IGNORE IT.
          IF(IBOUND(IC,IR,IL).LE.0) CYCLE
          C=EVTR(IC,IR)
          S=SURF(IC,IR)
          SS=S
          HH=HNEW(IC,IR,IL)
C
C--IF AQUIFER HEAD => SURF,SET Q=MAX ET RATE
          IF(HH.GE.SS) THEN
            QQ=-C
C
C--IF DEPTH=>EXTINCTION DEPTH, ET IS 0
C--OTHERWISE, LINEAR RANGE: Q=-HNEW*EVTR/EXDP -EVTR +EVTR*SURF/EXDP
          ELSE
            X=EXDP(IC,IR)
            XX=X
            DD=SS-HH
            IF(DD.GE.XX) THEN
              QQ=ZERO
            ELSE
              HHCOF=-C/X
              RRHS=(C*S/X)-C
              QQ= HH*HHCOF + RRHS
            ENDIF
          ENDIF
C
C--ADD Q TO BUFFER 1
          BUFF(IC,IR,1)=QQ
        ENDDO
      ENDDO
C
C--RECORD THEM.
      IF(NEVTOP.EQ.1) THEN
        IL=1
        IF(ILMTFMT.EQ.0) WRITE(IUMT3D)   ((IL,J=1,NCOL),I=1,NROW)
        IF(ILMTFMT.EQ.1) WRITE(IUMT3D,*) ((IL,J=1,NCOL),I=1,NROW)
      ELSEIF(NEVTOP.NE.1) THEN
        IF(ILMTFMT.EQ.0) WRITE(IUMT3D)   ((IEVT(J,I),J=1,NCOL),I=1,NROW)
        IF(ILMTFMT.EQ.1) WRITE(IUMT3D,*) ((IEVT(J,I),J=1,NCOL),I=1,NROW)
      ENDIF
C
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) ((BUFF(J,I,1),J=1,NCOL),I=1,NROW)
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) ((BUFF(J,I,1),J=1,NCOL),I=1,NROW)
      ENDIF
C
C--RETURN
      RETURN
      END
C
C
      SUBROUTINE LMT7GHB7(ILMTFMT,IUMT3D,KSTP,KPER,IGRID)
C *****************************************************************
C SAVE HEAD-DEPENDENT BOUNDARY CELL LOCATIONS AND VOLUMETRIC FLOW
C RATES FOR USE BY MT3D.
C *****************************************************************
C Modified from Harbaugh (2005)
C last modified: 08-08-2008
C
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,IBOUND,HNEW,BUFF
      USE GWFGHBMODULE,ONLY:NBOUND,BNDS
      CHARACTER*16 TEXT
      DOUBLE PRECISION CCGHB,CHB
C
C--SET POINTERS FOR THE CURRENT GRID     
c swm: already set in GWF2GHB7BD      CALL SGWF2GHB7PNT(IGRID)
C      
      TEXT='GHB'
      ZERO=0.      
C
C--WRITE AN IDENTIFYING HEADER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT,NBOUND
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT,NBOUND
      ENDIF
C
C--IF NO BOUNDARIES THEN SKIP
      IF(NBOUND.LE.0) GO TO 9999
C
C--FOR EACH GENERAL HEAD BOUND ACCUMULATE FLOW INTO AQUIFER
      DO L=1,NBOUND
C
C--GET LAYER, ROW AND COLUMN OF EACH GENERAL HEAD BOUNDARY.
        IL=BNDS(1,L)
        IR=BNDS(2,L)
        IC=BNDS(3,L)
C
C--RATE=0 IF IBOUND=<0
        RATE=ZERO
        IF(IBOUND(IC,IR,IL).GT.0) THEN
C
C--GET PARAMETERS FROM BOUNDARY LIST.          
          HB=BNDS(4,L)
          C=BNDS(5,L)
          CCGHB=C          
          CHB=C*HB
C
C--CALCULATE THE FOW RATE INTO THE CELL
          RATE= CHB - CCGHB*HNEW(IC,IR,IL)
        ENDIF
C
C--WRITE HEAD DEP. BOUND. LOCATION AND RATE
        IF(ILMTFMT.EQ.0) THEN
          WRITE(IUMT3D) IL,IR,IC,RATE
        ELSEIF(ILMTFMT.EQ.1) THEN
          WRITE(IUMT3D,*) IL,IR,IC,RATE
        ENDIF
      ENDDO
C
C--RETURN
 9999 RETURN
      END
C
C
      SUBROUTINE LMT7FHB7(ILMTFMT,IUMT3D,KSTP,KPER,IGRID)
C **********************************************************************
C SAVE SPECIFIED-FLOW CELL LOCATIONS AND VOLUMETRIC FLOW RATES
C FOR USE BY MT3D.
C **********************************************************************
C Modified from Leake and Lilly (1997), and Harbaugh (2005)
C last modified: 08-08-2008
C
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,IBOUND     
      USE GWFFHBMODULE,ONLY:NFLW,IFLLOC,BDFV
      CHARACTER*16 TEXT
C   
C--SET POINTERS FOR THE CURRENT GRID
c swm: already set in GWF2FHB7BD      CALL SGWF2FHB7PNT(IGRID)
C      
      TEXT='FHB'
      ZERO=0.
C
C--WRITE AN IDENTIFYING HEADER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT,NFLW
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT,NFLW
      ENDIF
C
C--IF NO SPECIFIED-FLOW CELL, RETURN
      IF(NFLW.LE.0) GO TO 9999
C
C--PROCESS SPECIFIED-FLOW CELLS ONE AT A TIME.

      DO L=1,NFLW
C
C--GET LAYER, ROW, AND COLUMN NUMBERS
        IR=IFLLOC(2,L)
        IC=IFLLOC(3,L)
        IL=IFLLOC(1,L)
        Q=ZERO
C
C--GET FLOW RATE FROM SPECIFIED-FLOW LIST
        IF(IBOUND(IC,IR,IL).GT.0) Q=BDFV(1,L)
C
C--WRITE SPECIFIED-FLOW CELL LOCATION AND RATE
        IF(ILMTFMT.EQ.0) THEN
          WRITE(IUMT3D) IL,IR,IC,Q
        ELSEIF(ILMTFMT.EQ.1) THEN
          WRITE(IUMT3D,*) IL,IR,IC,Q
        ENDIF
      ENDDO
C
C--NORMAL RETURN
 9999 RETURN
      END
C
C
      SUBROUTINE LMT7RES7(ILMTFMT,IUMT3D,KSTP,KPER,IGRID)
C **********************************************************************
C SAVE RESERVOIR CELL LOCATIONS AND VOLUMETRIC FLOW RATES
C FOR USE BY MT3D.
C **********************************************************************
C Modified from Fenske et al., (1996), Harbaugh (2005)
C last modified: 08-08-2008
C
      USE GLOBAL ,      ONLY: HNEW,IBOUND,BUFF,NCOL,NROW,NLAY 
      USE GWFRESMODULE, ONLY: NRES,NRESOP,IRES,IRESL,BRES,CRES,
     &                        BBRES,HRES  
      CHARACTER*16 TEXT    
C
C--SET POINTERS FOR THE CURRENT GRID      
c swm: already set in GWF2RES7BD      CALL SGWF2RES7PNT(IGRID)            
C      
      TEXT='RES'
      ZERO=0.
C
C--CLEAR BUFFER
      DO IL=1,NLAY
        DO IR=1,NROW
          DO IC=1,NCOL
            BUFF(IC,IR,IL)=ZERO
          ENDDO
        ENDDO
      ENDDO
C
C--FOR EACH RESERVOIR REACH ACCUMULATE RESERVOIR FLOW
      DO 200 I=1,NROW
      DO 190 J=1,NCOL
      NR=IRES(J,I)
      IF(NR.LE.0) GO TO 190
      IF(NR.GT.NRES) GO TO 190
      IR=I
      IC=J
C
C--FIND LAYER NUMBER FOR RESERVOIR CELL
      IF(NRESOP.EQ.1) THEN
       IL=1
      ELSE IF(NRESOP.EQ.2) THEN
       IL=IRESL(IC,IR)
      ELSE
       DO 60 K=1,NLAY
       IL=K
C--UPPERMOST ACTIVE CELL FOUND, SAVE LAYER INDEX IN 'IL'
       IF(IBOUND(IC,IR,IL).GT.0) GO TO 70
C--SKIP THIS CELL IF VERTICAL COLUMN CONTAINS A CONSTANT-
C--HEAD CELL ABOVE RESERVOIR LOCATION
       IF(IBOUND(IC,IR,IL).LT.0) GO TO 190
   60  CONTINUE
       GO TO 190
      ENDIF
C
C--IF THE CELL IS EXTERNAL SKIP IT.
      IF(IBOUND(IC,IR,IL).LE.0) GO TO 190
C
C--IF RESERVOIR STAGE IS BELOW RESERVOIR BOTTOM, SKIP IT
   70 HR=HRES(NR)
      IF(HR.LE.BRES(IC,IR))  GO TO 190
C--SINCE RESERVOIR IS ACTIVE AT THIS LOCATION,
C--GET THE RESERVOIR DATA.
      CR=CRES(IC,IR)
      RBOT=BBRES(IC,IR)
      HHNEW=HNEW(IC,IR,IL)
C
C--COMPUTE RATE OF FLOW BETWEEN GROUND-WATER SYSTEM AND RESERVOIR.
C
C--GROUND-WATER HEAD > BOTTOM THEN RATE=CR*(HR-HNEW).
      IF(HHNEW.GT.RBOT) RATE=CR*(HR-HHNEW)
C
C--GROUND-WATER HEAD < BOTTOM THEN RATE=CR*(HR-RBOT)
      IF(HHNEW.LE.RBOT) RATE=CR*(HR-RBOT)
C
C--ADD RATE TO BUFFER.
      BUFF(IC,IR,IL)=BUFF(IC,IR,IL)+RATE
  190 CONTINUE
  200 CONTINUE
C
C--COUNT RES CELLS WITH NONZERO FLOW RATE
      NTEMP=0
      DO IL=1,NLAY
        DO IR=1,NROW
          DO IC=1,NCOL
            IF(IBOUND(IC,IR,IL).LE.0) CYCLE
            IF(BUFF(IC,IR,IL).NE.ZERO) THEN
              NTEMP=NTEMP+1
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--WRITE AN IDENTIFYING HEADER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT,NTEMP
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT,NTEMP
      ENDIF
C
C--IF NO RES CELLS WITH NONZERO Q, RETURN
      IF(NTEMP.EQ.0) GO TO 9999
C
C--WRITE RES CELL LOCATION AND FLOW RATE
      DO IL=1,NLAY
        DO IR=1,NROW
          DO IC=1,NCOL
            IF(IBOUND(IC,IR,IL).LE.0) CYCLE
            RATE=BUFF(IC,IR,IL)
            IF(RATE.NE.ZERO) THEN
              IF(ILMTFMT.EQ.0) THEN
                WRITE(IUMT3D)   IL,IR,IC,RATE
              ELSEIF(ILMTFMT.EQ.1) THEN
                WRITE(IUMT3D,*) IL,IR,IC,RATE
              ENDIF  
            ENDIF
          ENDDO
        ENDDO
      ENDDO
C
C--NORMAL RETURN
 9999 RETURN
      END
C
C
      SUBROUTINE LMT7STR7(ILMTFMT,IUMT3D,KSTP,KPER,IGRID)
C **********************************************************************
C SAVE STREAM CELL LOCATIONS AND VOLUMETRIC FLOW RATES FOR USE BY MT3D.
C **********************************************************************
C Modified from Prudic (1989), Harbaugh (2005)
C last modified: 08-08-2008
C
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,IBOUND      
      USE GWFSTRMODULE,ONLY:NSTREM,STRM,ISTRM
      CHARACTER*16 TEXT     
C
C--SET POINTERS FOR THE CURRENT GRID      
c swm: already set in GWF2STR7BD   CALL SGWF2STR7PNT(IGRID)    
C                    
      TEXT='STR'
      ZERO=0.
C
C--WRITE AN IDENTIFYING HEADER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT,NSTREM
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT,NSTREM
      ENDIF
C
C--IF NO REACHES, SKIP
      IF(NSTREM.EQ.0) GO TO 9999
C
C--FOR EACH STREAM REACH GET LEAKAGE TO OR FROM IT
      DO L=1,NSTREM
C
C--GET REACH LOCATION AND FLOW RATE
        IL=ISTRM(1,L)
        IR=ISTRM(2,L)
        IC=ISTRM(3,L)
        IF(IBOUND(IC,IR,IL).LE.0) THEN
          RATE=ZERO
        ELSE
          RATE=STRM(11,L)
        ENDIF
C
C--WRITE STREAM REACH LOCATION AND RATE
        IF(ILMTFMT.EQ.0) THEN
          WRITE(IUMT3D) IL,IR,IC,RATE
        ELSEIF(ILMTFMT.EQ.1) THEN
          WRITE(IUMT3D,*) IL,IR,IC,RATE
        ENDIF
C        
      ENDDO
C
C--NORMAL RETURN
 9999 RETURN
      END
C
C
      SUBROUTINE LMT7MNW27(ILMTFMT,IUMT3D,KSTP,KPER,IGRID)
C *********************************************************************
C SAVE MNW LOCATIONS AND VOLUMETRIC FLOW RATES FOR USE BY MT3D.
C *********************************************************************
C Modified from MNW by Halford and Hanson (2002)
C last modification: 08-08-2008
C Modified from MNW2 by Konikow and Hornberger (2009)
C modification: 10-21-2010:  swm  
C last modification: 2-16-2012:  awh
C
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,IBOUND
      USE GWFMNW2MODULE,ONLY:NMNW2,NTOTNOD,MNW2,MNWNOD,MNWMAX
      INTEGER firstnode, lastnode
      CHARACTER*16 TEXT
C
C--SET POINTERS FOR THE CURRENT GRID
c swm: already set in GWF2MNW7BD      CALL SGWF2MNW7PNT(IGRID)
C      
      TEXT='MNW'
      ZERO=0.
c swm: SET NUMBER OF ACTIVE WELL NODES BASED ON NMNW2 AND NTOTNOD
      NACTW=NTOTNOD
      IF(NMNW2.LE.0) NACTW=0
C
C--WRITE AN IDENTIFYING HEADER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT,NACTW
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT,NACTW
      ENDIF
C
C--IF THERE ARE NO WELLS RETURN
      IF(NMNW2.LE.0) GO TO 9999
C
C--PROCESS WELL LIST -- write Q for all nodes even if inactive
c  Loop over all wells
      DO iw=1,MNWMAX
        IDwell = iw
c  active well check
        firstnode=MNW2(4,iw)
        lastnode=MNW2(4,iw)+ABS(MNW2(2,iw))-1
c   Loop over nodes in well
        do INODE=firstnode,lastnode
          il=MNWNOD(1,INODE)              
          ir=MNWNOD(2,INODE)              
          ic=MNWNOD(3,INODE)              
C
C--IF CELL IS EXTERNAL OR WELL INACTIVE Q=0
          IF(IBOUND(IC,IR,IL).LE.0 .OR. MNW2(1,iw).LE.0.) then
            Q=ZERO
          else
            Q=MNWNOD(4,INODE)
          end if
C
C--DUMMY VARIABLE QSW NOT USED, SET TO 0
          QSW=ZERO
C
C--SAVE TO OUTPUT FILE
          IF(ILMTFMT.EQ.0) THEN
            WRITE(IUMT3D) IL,IR,IC,Q,IDwell,QSW
          ELSEIF(ILMTFMT.EQ.1) THEN
            WRITE(IUMT3D,*) IL,IR,IC,Q,IDwell,QSW
          ENDIF
        enddo
      ENDDO
C
C--RETURN
 9999 RETURN
      END
C
      SUBROUTINE LMT7MNW17(ILMTFMT,IUMT3D,KSTP,KPER,IGRID)
C *********************************************************************
C SAVE MNW LOCATIONS AND VOLUMETRIC FLOW RATES FOR USE BY MT3D.
C *********************************************************************
C Modified from MNW by Halford and Hanson (2002)
C last modification: 08-08-2008
C swm: 02-12-2010
C
      USE GLOBAL,      ONLY:NCOL,NROW,NLAY,IBOUND
      USE GWFMNW1MODULE,ONLY:NWELL2,WELL2
      CHARACTER*16 TEXT
C
C--SET POINTERS FOR THE CURRENT GRID
cswm:      CALL SGWF2MNW1PNT(IGRID)
C      
      TEXT='MNW'
      ZERO=0.
C
C--WRITE AN IDENTIFYING HEADER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT,NWELL2
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT,NWELL2
      ENDIF
C
C--IF THERE ARE NO WELLS RETURN
      IF(NWELL2.LE.0) GO TO 9999
C
C--PROCESS WELL LIST
      DO m = 1,nwell2
        n = ifrl( well2(1,m) )
        il = (n-1) / (ncol*nrow) + 1
        ir = mod((n-1),ncol*nrow)/ncol + 1
        ic = mod((n-1),ncol) + 1
        IDwell = ifrl(well2(18,m))  !IDwell in well2(18,m); cdl 4/19/05
        Q = well2(17,m)
C
C--IF CELL IS EXTERNAL Q=0
        IF(IBOUND(IC,IR,IL).LE.0) Q=ZERO
C
C--DUMMY VARIABLE QSW NOT USED, SET TO 0
        QSW=ZERO
C
C--SAVE TO OUTPUT FILE
        IF(ILMTFMT.EQ.0) THEN
          WRITE(IUMT3D) IL,IR,IC,Q,IDwell,QSW
        ELSEIF(ILMTFMT.EQ.1) THEN
          WRITE(IUMT3D,*) IL,IR,IC,Q,IDwell,QSW
        ENDIF
      ENDDO
C
C--RETURN
 9999 RETURN
      END
C
      SUBROUTINE LMT7ETS7(ILMTFMT,IUMT3D,KSTP,KPER,IGRID)
C ********************************************************************
C SAVE SEGMENTED EVAPOTRANSPIRATION LAYER INDICES (IF NLAY>1) AND
C VOLUMETRIC FLOW RATES FOR USE BY MT3D.
C ********************************************************************
C Modified from Banta (2000), Harbaugh (2005)
C last modified: 08-08-2008
C
      USE GLOBAL ,      ONLY: HNEW,IBOUND,BUFF,NCOL,NROW,NLAY      
      USE GWFETSMODULE, ONLY: NETSOP,NETSEG,IETS,ETSR,ETSX,ETSS,
     &                        PXDP,PETM
      CHARACTER*16 TEXT
      DOUBLE PRECISION QQ,HH,SS,DD,XX,HHCOF,RRHS,PXDP1,PXDP2
C
C--SET POINTERS FOR THE CURRENT GRID
c swm: already set in GWF2ETS7BD      CALL SGWF2ETS7PNT(IGRID)
C      
      TEXT='ETS'
      ZERO=0.      
C
C--WRITE AN IDENTIFYING HEADER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT
      ENDIF      
C
C--CLEAR THE BUFFER
      DO IL=1,NLAY
        DO IR=1,NROW
          DO IC=1,NCOL
            BUFF(IC,IR,IL)=ZERO
          ENDDO   
        ENDDO   
      ENDDO   
C
C--PROCESS EACH HORIZONTAL CELL LOCATION
      DO IR=1,NROW
        DO IC=1,NCOL
C
C--SET THE LAYER INDEX EQUAL TO 1.
          IL=1
C
C--IF OPTION 2 IS SPECIFIED THEN GET LAYER INDEX FROM IETS ARRAY
          IF (NETSOP.EQ.2) IL=IETS(IC,IR)
          IF (IL.EQ.0) CYCLE
C
C--IF CELL IS EXTERNAL THEN IGNORE IT.
          IF (IBOUND(IC,IR,IL).LE.0) CYCLE
C          
          C=ETSR(IC,IR)
          S=ETSS(IC,IR)
          SS=S
          HH=HNEW(IC,IR,IL)
C
C--IF HEAD IN CELL => ETSS,SET Q=MAX ET RATE.
          IF (HH.GE.SS) THEN
            QQ=-C
          ELSE
C
C--IF DEPTH=>EXTINCTION DEPTH, ET IS 0.
            X=ETSX(IC,IR)
            XX=X
            DD=SS-HH
            IF (DD.LT.XX) THEN
C--VARIABLE RANGE.  CALCULATE Q DEPENDING ON NUMBER OF SEGMENTS
C
              IF (NETSEG.GT.1) THEN
C               DETERMINE WHICH SEGMENT APPLIES BASED ON HEAD, AND
C               CALCULATE TERMS TO ADD TO RHS AND HCOF
C
C               SET PROPORTIONS CORRESPONDING TO ETSS ELEVATION
                PXDP1 = 0.0
                PETM1 = 1.0
                DO ISEG = 1,NETSEG
C                 SET PROPORTIONS CORRESPONDING TO LOWER END OF
C                 SEGMENT
                  IF (ISEG.LT.NETSEG) THEN
                    PXDP2 = PXDP(IC,IR,ISEG)
                    PETM2 = PETM(IC,IR,ISEG)
                  ELSE
                    PXDP2 = 1.0
                    PETM2 = 0.0
                  ENDIF
                  IF (DD.LE.PXDP2*XX) THEN
C                   HEAD IS IN DOMAIN OF THIS SEGMENT
                    EXIT
                  ENDIF
C                 PROPORTIONS AT LOWER END OF SEGMENT WILL BE FOR
C                 UPPER END OF SEGMENT NEXT TIME THROUGH LOOP
                  PXDP1 = PXDP2
                  PETM1 = PETM2
                ENDDO   
C--CALCULATE ET RATE BASED ON SEGMENT THAT APPLIES AT HEAD
C--ELEVATION
                HHCOF = -(PETM1-PETM2)*C/((PXDP2-PXDP1)*X)
                RRHS = -HHCOF*(S-PXDP1*X) - PETM1*C
              ELSE
C--SIMPLE LINEAR RELATION.  Q=-ETSR*(HNEW-(ETSS-ETSX))/ETSX, WHICH
C--IS FORMULATED AS Q= -HNEW*ETSR/ETSX + (ETSR*ETSS/ETSX -ETSR).
                HHCOF = -C/X
                RRHS = (C*S/X) - C
              ENDIF
              QQ = HH*HHCOF + RRHS
            ELSE
              QQ = ZERO
            ENDIF
          ENDIF  
C
C--ADD Q TO BUFFER.
          Q=QQ
          BUFF(IC,IR,1)=Q
        ENDDO   
      ENDDO   
C
C--RECORD THEM
      IF(NETSOP.EQ.1) THEN
        IL=1
        IF(ILMTFMT.EQ.0) WRITE(IUMT3D)  ((IL,J=1,NCOL),I=1,NROW)
        IF(ILMTFMT.EQ.1) WRITE(IUMT3D,*)((IL,J=1,NCOL),I=1,NROW)
      ELSEIF(NETSOP.NE.1) THEN
        IF(ILMTFMT.EQ.0) THEN 
          WRITE(IUMT3D)   ((IETS(J,I),J=1,NCOL),I=1,NROW)
        ELSEIF(ILMTFMT.EQ.1) THEN
          WRITE(IUMT3D,*) ((IETS(J,I),J=1,NCOL),I=1,NROW)
        ENDIF  
      ENDIF
C
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D)   ((BUFF(J,I,1),J=1,NCOL),I=1,NROW)
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) ((BUFF(J,I,1),J=1,NCOL),I=1,NROW)
      ENDIF
C
C--RETURN
      RETURN
      END
C
C
      SUBROUTINE LMT7DRT7(ILMTFMT,IUMT3D,KSTP,KPER,IGRID)
C ******************************************************************
C SAVE DRT (Drain with Return Flow) CELL LOCATIONS AND 
C VOLUMETRIC FLOW RATES FOR USE BY MT3D.
C ******************************************************************
C Modified from Banta (2000), Harbaugh (2005)
C last modified: 08-08-2008
C
      USE GLOBAL ,      ONLY: HNEW,IBOUND,NCOL,NROW,NLAY
      USE GWFDRTMODULE, ONLY: DRTF,NDRTCL,IDRTFL,NRFLOW
      CHARACTER*16 TEXT
      DOUBLE PRECISION HHNEW,EEL,CC,CEL,QQ,QQIN
C
C--SET POINTERS FOR THE CURRENT GRID
c swm: already set in GWF2DRT7BD      CALL SGWF2DRT7PNT(IGRID)
C      
      TEXT='DRT'
      ZERO=0.
C
C--WRITE AN IDENTIFYING HEADER
      IF(ILMTFMT.EQ.0) THEN
        WRITE(IUMT3D) KPER,KSTP,NCOL,NROW,NLAY,TEXT,NDRTCL+NRFLOW
      ELSEIF(ILMTFMT.EQ.1) THEN
        WRITE(IUMT3D,*) KPER,KSTP,NCOL,NROW,NLAY
        WRITE(IUMT3D,*) TEXT,NDRTCL+NRFLOW
      ENDIF      
C
C--IF THERE ARE NO DRAIN-RETURN CELLS, SKIP.
      IF (NDRTCL+NRFLOW.LE.0) GO TO 9999
C
C--LOOP THROUGH EACH DRAIN-RETURN CELL, CALCULATING FLOW.
      DO L=1,NDRTCL
C
C--GET LAYER, ROW & COLUMN OF CELL CONTAINING DRAIN.
        IL=DRTF(1,L)
        IR=DRTF(2,L)
        IC=DRTF(3,L)
        Q=ZERO
        ILR=0
        IF(IDRTFL.GT.0) THEN
          QIN=ZERO
          ILR=DRTF(6,L)
          IRR=DRTF(7,L)
          ICR=DRTF(8,L)
          IF(IBOUND(ICR,IRR,ILR).LE.0) ILR=0
        ENDIF                
C
C--IF CELL IS NO-FLOW OR CONSTANT-HEAD, IGNORE IT.
        IF (IBOUND(IC,IR,IL).LE.0) GOTO 99
C
C--GET DRAIN PARAMETERS FROM DRAIN-RETURN LIST.
        EL=DRTF(4,L)
        EEL=EL
        C=DRTF(5,L)
        HHNEW=HNEW(IC,IR,IL)
C
C--IF HEAD HIGHER THAN DRAIN, CALCULATE Q=C*(EL-HHNEW).
        IF(HHNEW.GT.EEL) THEN
          CC=C
          CEL=C*EL
          QQ=CEL - CC*HHNEW
          Q=QQ          
          IF(IDRTFL.GT.0) THEN           
            IF(ILR.NE.0) THEN
              RFPROP = DRTF(9,L)
              QQIN = RFPROP*(CC*HHNEW-CEL)
              QIN = QQIN
            ENDIF
          ENDIF
        ENDIF
   99   CONTINUE     
C
C--WRITE DRT LOCATION AND RATE (both host and recipient)
        mhost=0
        QSW=ZERO
C       main drain (host to recipient cell)
        IF(ILMTFMT.EQ.0) WRITE(IUMT3D)   IL,IR,IC,Q,mhost,QSW
        IF(ILMTFMT.EQ.1) WRITE(IUMT3D,*) IL,IR,IC,Q,mhost,QSW 
C       return flow recipient cell 
        if(IDRTFL.GT.0 .AND. ILR.GT.0) then
          mhost = ncol*nrow*(IL-1) + ncol*(IR-1) + IC
          IF(ILMTFMT.EQ.0) THEN
            WRITE(IUMT3D)   ILR,IRR,ICR,QIN,mhost,QSW            
          ELSEIF(ILMTFMT.EQ.1) THEN 
            WRITE(IUMT3D,*) ILR,IRR,ICR,QIN,mhost,QSW
          ENDIF
        endif
      ENDDO   
C
C--RETURN
 9999 RETURN
      END
C***********************************************************************
      SUBROUTINE SLMT7PNT(IGRID)
C     ******************************************************************
C     CHANGE POINTERS FOR LMT DATA TO A DIFFERENT GRID
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE LMTMODULE
C     ------------------------------------------------------------------
C
      ISSMT3D=>LMTDAT(IGRID)%ISSMT3D
      IUMT3D=>LMTDAT(IGRID)%IUMT3D
      ILMTFMT=>LMTDAT(IGRID)%ILMTFMT
      RETURN
      END
C***********************************************************************
      SUBROUTINE SLMT7PSV(IGRID)
C     ******************************************************************
C     SAVE POINTERS ARRAYS FOR LMT DATA TO APPROPRIATE GRID
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE LMTMODULE
C     ------------------------------------------------------------------
C
      LMTDAT(IGRID)%ISSMT3D=>ISSMT3D
      LMTDAT(IGRID)%IUMT3D=>IUMT3D
      LMTDAT(IGRID)%ILMTFMT=>ILMTFMT
C
      RETURN
      END
C***********************************************************************
      SUBROUTINE LMT7DA(IGRID)
C     ******************************************************************
C     DEALLOCATE LMT DATA
C     ******************************************************************
C
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE LMTMODULE
C     ------------------------------------------------------------------
C
      DEALLOCATE(LMTDAT(IGRID)%ISSMT3D)
      DEALLOCATE(LMTDAT(IGRID)%IUMT3D)
      DEALLOCATE(LMTDAT(IGRID)%ILMTFMT)
C
      RETURN
      END
